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Abstract 

(N 

,_^ We propose a novel compressed sensing technique to accelerate the magnetic resonance imaging (MRI) acquisition 

^^ process. The method, coined spread spectrum MRI or simply S2MRI, consists of pre-modulating the signal of interest 

by a linear chirp before random fc-space under-sampling, and then reconstructing the signal with non-linear algorithms 
^ Cy that promote sparsity. The effectiveness of the procedure is theoretically underpinned by the optimization of the 

^ .^ coherence between the sparsity and sensing bases. The proposed technique is thoroughly studied by means of numerical 

^\ simulations, as well as phantom and in vivo experiments on a 7T scanner. Our results suggest that S2MRI performs 

better than state-of-the-art variable density fe-space under-sampling approaches. 

X 

^^ Index Terms 

(^ compressed sensing, spread spectrum, magnetic resonance imaging. 

^ I. Introduction 

in 

C^ In magnetic resonance imaging (MRI), the signal of interest p represents the magnetization induced by resonance 

p^ in the imaged tissues. In the standard setting, data acquired in MRI provide complete Fourier, or fc-space, measure- 

ff") ments of this signal p. Accelerating the acquisition process, or equivalently increasing the achievable resolution for 

o 

p,i a fixed acquisition time, is of major interest for MRI applications. Recent approaches based on compressed sensing 

. . seek to reconstruct the signal from incomplete fc-space information, hence defining an ill-posed inverse problem. 

> 

• ^H The ill-posed inverse problem is regularized by the introduction of sparsity priors, acknowledging the fact that many 

}_j MRI signals are sparse in well-chosen bases; i.e., that their expansion contains only a small number of non-zero 

coefficients. Such compressed sensing techniques have been developed for static and dynamic imaging [l]-[5]. 
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parallel MRI [6]-[8], MR spectroscopic imaging [9]-[ll], and many other applications. Several algorithms have 
also been proposed to reconstruct MRI signals from under-sampled fc-space measurements (see, e.g., [12]-[14]). 

In the framework of compressed sensing, signals are usually measured through random matrices to ensure that 
any sparse signal can be recovered with overwhelming probability. The common approach in MRI consists in also 
exploiting the fact that the energy of MRI signals is usually concentrated at low spatial frequencies. Therefore a 
variable density fc-space random sampling where the under-sampling ratio increases at high frequencies is usually 
used [1]. This method, heuristic in nature, was shown to be very effective in enhancing the signal reconstruction 
quality when random distributions optimizing the associated point spread function are used. Other approaches that 
optimize the acquisition procedure have also been proposed: fc-space sampling optimization by Bayesian inference 
[15]; random fc-space convolution with Toeplitz matrices [16]-[18]; encoding by projection onto random waveforms 
with Gaussian distributions [19]. 

In the present work, we propose the use of a spread spectrum technique to accelerate single coil MRI acquisition 
in the framework of compressed sensing. We study this method, coined spread spectrum MRI or simply S2MRI, 
theoretically, numerically via simulations, and empirically via real acquisitions. The essence of our strategy consists 
of pre-modulating the image by a linear chirp, which results from the application of quadratic phase profiles, and then 
performing random fc-space under-sampling. Images are then reconstructed with non-linear algorithms promoting 
signal sparsity. The enhancement of the signal reconstruction quality is linked to a decrease of coherence of the 
measurement system [20], [21]. In MRI, this type of modulation is known as phase scrambling. It can be obtained 
by using dedicated coils or by modifying radio frequency (RF) pulses. It has been used for various purposes, such 
as improving dynamic range [22], [23], or reducing aliasing artifacts [24], [25], but never in a compressed sensing 
perspective. 

Let us acknowledge that this spread spectrum technique was initially introduced by some of the authors for 
compressive sampling of pulse trains in [26]. Its transfer to a setting encompassing analog signals and modulations 
was studied in the context of radio interferometry [21], [27], [28]. The effectiveness of the method for MRI was 
briefly discussed in [29]-[31]. 

The paper is organized as follows. In Section II, we explain the principle of the spread spectrum technique 
in a simplified analog setting. In Section III, we formulate the inverse problem for image reconstruction from 
under-sampled fc-space measurements in the presence of chirp modulation and compare the S2MRI technique to the 
variable density sampling method on the basis of numerical simulations of MR acquisitions of a brain image. In 
Section IV, we describe our implementation of the chirp modulation (on a 7T scanner) and show the effectiveness 
of S2MRI using numerical simulations in this precise setting and real acquisitions of phantom and in vivo data. 
Finally, we conclude and discuss potential evolutions of the technique in Section V. 
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II. Spread spectrum principle 

A. Compressed sensing 

The essence of the recent theory of compressed sensing is the merging of data acquisition and compression [32]- 
[39]. Beyond MRI, it is well-known that a large variety of natural signals are sparse or compressible in multi-scale 
bases, such as wavelet bases. In the compressed sensing theory, signals are usually expressed as A^-dimensional 
vectors: p e C^. By definition, a signal is sparse in some orthonormal basis U^ e C^^^, called a sparsity basis, 
if its expansion contains only a small number K -^ N of non-zero coefficients. More generally it is compressible 
if its expansion contains only a small number of significant coefficients. The decomposition of p in U^ is denoted 
a e C^, and it satisfies 

p^^a. (1) 

The theory of compressed sensing demonstrates that a small number M ^ iV of linear and non-adaptative 
measurements u E C^^ suffices for an accurate and stable reconstruction of the signal p. These measurements may, 
for example, be obtained by projection onto M randomly selected basis vectors of an orthonormal basis <t> G C^^^, 
called a sensing basis. This selection process can be modeled by a multiplication with a rectangular binary matrix 
M G M*^><^ that contains only one non-zero value on each line, at the index of the basis vector to be selected. 
To model a non-perfect sensing process, these measurements are assumed to be contaminated by independent and 
identically distributed noise n E C^^. The measurement model thus satisfies 

v> = ea + n, with 9 = M0*v|/ e M^^^^, (2) 

where the matrix identifies the measurement matrix as seen from the sparsity basis (•* denotes the conjugate 
transpose operation). 

To reconstruct the signal p from the measurements v, the compressed sensing framework proposes, among 
other approaches, to solve the Basis Pursuit (BP) minimization problem [32]-[39]. This problem regularizes the 
originally ill-posed inverse problem related to (2) with an explicit sparsity or compressibility prior on the signal. 
In the presence of noise, the BP problem is the minimization of the £i norm' of a under a constraint on the £2 
norm^ of the residual noise fi ~ v Qa: 

a* = arg min ||q:||i subject to \\v — 0ci||2 < £• (3) 

The corresponding reconstructed signal is p* = '^a*. 

In this setting, the compressed sensing theory shows that if the number of measurements satisfies 

M > DNK^? (<t>* , M/) In'' TV, (4) 

for a universal constant D, then a. is recovered accurately by solving the BP problem (3) [39]. In relation (4), 
ji (<J>*, U') is the mutual coherence between the sparsity and sensing bases. It is defined as the maximum projection, 

'||q:||i = J^j_i \cti\ where ai is the ith entry of the vector a. 
^||n||2 = (Yl^—i |"iP)'''^ where n^ is the ith entry of the vector n. 
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Fig. 1. Spread spectrum principle. A signal p (top left panel) is modulated by a linear chirp (top right panel). In fc-space, the modulation 
amounts to the convolution of the Fourier transform of the signal (middle left panel) with that of the chirp (middle right panel). The spectrum 
of the resulting signal (bottom panel) spreads out in fe-space. 



in absolute value, between the sparsity basis vectors (pi, I < i < N, and sensing bases vectors xpi', 1 < i' < N 
[36], [39]: 

iV-i/2<^(0*,i|/)= max \<f^-ip^'\<l. (5) 

l<i,i'<N 

The mutual coherence fj,{'t>*,'^) plays a crucial role in relation (4). Indeed, the number of measurements M* 
needed to reconstruct X-sparse signals increases quadratically with its value: M* ex /i^(0*, U'). In the worst case 
where /i (0*, U') = 1, M* is of the order of the signal dimension N and under-sampling is impossible. However, 
when the mutual coherence is at its minimum, N~^^'^, M* is reduced to the order of the sparsity level K. This result 
can intuitively be explained with a simple consideration on the spread of the energy of the sparsity basis vectors in 
the measurement domain 0. In an incoherent orthonormal system, the absolute value of the scalar product between 
the sensing basis vector ipt and the sparsity basis vector ■0^/ is small for all pairs of indices 1 < i,i' < A^. As 
II 't'*''/'?:' II 2 = ll^'/'j'lb- the energy of the sparsity basis vector ■0^/ spreads equally over the sensing basis vectors tp^. 
Consequently, whatever index i is selected to perform a measurement, one always gets information concerning all 
the sparsity basis vectors describing the original signal. The number of measurements needed for accurate recovery 
thus decreases. 
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TABLE I 

Values of A^c ^il at different chirp rates w 



N.t^l 


Dirac 


Haar 


Fourier 


tD = 


1.00 


256 


256 


to = 0.1 


2.27 


43.5 


15.5 


to = 0.3 


2.97 


25.9 


6.11 


tD = 0.5 


3.29 


22.4 


4.17 



B. Spread spectrum and coherence reduction 

The above considerations suggest that to reduce to the number of measurements needed for accurate recovery of 
MRI signals, one can try to modify the MRI acquisition procedure in such a way that the energy of the sparsity basis 
vectors spreads all over the fc-space. Following this idea, the S2MRI technique introduces a linear chirp modulation 
of the signal of interest before random selection of fc-space coefficients (see Fig. 1). Indeed, this modulation 
conserves the energy of the input signal and corresponds to a convolution that generically spreads the spectrum. 

To study theoretically the proposed acquisition scheme, we consider a simplified analog setting where the one- 
dimensional signal of interest is denoted by a complex-valued function p{x) of the position a; e M. This signal is 
limited on a finite field of view L. For the sake of simplicity of the following theoretical result only, we assume 
that the energy of this signal beyond the spatial frequency B is negligible^. We thus sample this signal on a discrete 
uniform grid of iV = 2LB points and represent the resulting discrete signal by a vector p E C^. The vector p is 
assumed to be ii'-sparse in an orthonormal basis U' e C^^^. We also consider a one-dimensional Unear chirp, with 
chirp rate w G M, that reads as a complex-valued function c{x) = 0™"*^ . On the field of view L, this linear chirp is 
approximately band limited at its maximum instantaneous frequency \w\L/2. This band limit can be parametrized 
in terms of a discrete chirp rate w — wL'^ /N and thus \w\L/2 ~ \w\B. 

In this setting, the S2MRI measurement model is given by equation (2) with 

0* = F*CU e C^-^^^. (6) 

In the above equation, the matrix U represents an up-sampling operator needed to avoid aliasing of the modulated 
signal due to a lack of sampling resolution in a digital description of the originally analog problem. Indeed, the 
convolution in Fourier space induced by the analog chirp modulation implies that the band limit of the modulated 
signal is the sum of the individual band limits of the original signal and of the chirp c. Therefore, an up-sampled 
grid with at least N^. = {l + \w\)N points needs to be considered and the modulated signal is correctly obtained by 
applying the chirp modulation on the signal after up-sampling on the Nc points grid. The up-sampling operator U, 
implemented in Fourier space by zero padding, is thus of size C^"^'^ and satisfies U*U = I G C^^^ . The matrix 
C e C^cXAf,, jg jjjg diagonal matrix implementing the chirp modulation on the up-sampled grid and the matrix 

^In section III, we will introduce a setting that properly takes into account the fact that, in general, MRI signals are not band limited. 
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Fig. 2. Probability of recovery of the signal p as a function of the number of measurements M for to = (dotted black curve), w = 0.1 
(continuous blue curve), to = 0.3 (continuous red curve), w = 0.5 (continuous black curve), and three sparsity bases: the Dirac basis (left); 
the Haar wavelet basis (middle); the Fourier basis (right). 



F = {fi\\<i<is! Stands for the discrete Fourier basis on the same grid. 

The S2MRI sensing matrix is not orthogonal because of the presence of the matrix U. Consequently, the recovery 
condition (4) does not strictly hold. However, one can obtain a similar recovery condition [20]. In particular, if 



M > DNcfii{<i>*,^')K\og'{N), 



(7) 



for some universal constant D, then the vector p is accurately recovered with high probability by solving the BP 
problem (3). This result shows that the number of measurements M needed to reconstruct A'-sparse signals is 
proportional to the product iVc /i^(0*, U') rather than to /i^(0*,U') only. The S2MRI technique is efficient only 
if this product decreases with the chirp rate w. In other words, the number of measurements needed for accurate 
recovery of sparse signals decreases only for sparsity basis U^ for which the mutual coherence /x^(<J>*, U^) decreases 
faster than N^ increases. 

C. Illustration 

We illustrate here the effect of the chirp modulation on the number of measurements needed for accurate recovery 
of If -sparse signals. 

We consider a 1-dimensional complex signal p of size N — 256 corresponding to one line of an MRI brain 
image. This signal is decomposed into three different sparsity bases U' and hard-thresholded at a fixed sparsity 
K — 25. The sparsity bases considered are the Dirac basis, the Haar wavelet basis and the Fourier basis. This 
signal is then probed according to relation (2) with 0* — F*CU and w G {0, 0.1, 0.3, 0.5}, and reconstructed from 
different numbers of measurements M by solving the BP problem (3). Each time, the probability of recovery "* of 
the signal is computed over 1000 simulations. In this experiment, no noise is added to the measurements y, and 
the indices i of the selected Fourier basis vectors fi are chosen uniformly at random from {0, . . . , A^c ^ 1}- 

The BP reconstructions in the presence of chirp modulation with chirp rate w are denoted by BPw. For each 
sparsity basis and chirp rate considered, the curves of the probability of recovery as a function of the number of 

'^Perfect recovery is considered to occur if the £2-norm between the original signal p and the reconstructed signal p* satisfies: ||p — p*||2 < 
10-^i|p||2 
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measurements M are reported in Fig. 2. The corresponding values of the product N^nli'i'* ,'^) are reported in 
Table I. 

For the Dirac basis, the number of measurements needed to reach a probability of recovery of 1 slightly increases 
with the chirp rate, as suggested by the values of iVc /i^(4>*, U^) in Table I and relation (7). On the contrary, for 
the Haar and Fourier bases, the values of A^c M?!*^'*; ^) ir* the table predict a drastic improvement of the results 
with an increase of the chirp rate. This prediction is confirmed by the results presented in Fig. 2. Note that for the 
Haar and Fourier bases, the value of the product N^ iil('i>* , U') decreases much less between w — 0.3 and w — 0.5 
than between w ~ 0.1 and w ~ 0.3, suggesting a smaller improvement in the number of measurements needed 
for accurate recovery. This is also in line with the curves of the probability of recovery in Fig. 2. In summary, as 
predicted by the theory, the effect of the S2MRI technique depends on the sparsity basis. Note that the decrease of 
the performance for the Dirac basis (optimally incoherent only at w) = 0) is negligible compared to the improvement 
obtained for the two other bases. Note also that MRI signals are usually sparse in wavelet bases [1]. The results 
obtained with the Haar wavelet basis therefore suggest strong efficiency of the technique in MRI. 

One can wonder if the proposed encoding scheme can result in shift-variant reconstruction quality. Indeed, the 
phase variation at the center of the chirp is less rapid than at the limit of the field of view. Let us consider the 
case of a signal sparse in the Haar wavelet basis. This signal can, roughly, be separated in large scale and fine 
scale sparsity basis vectors. The region where the chirp is oscillating slowly is limited to a small part of the field 
of view. Consequently, large scale sparsity basis vectors are necessarily affected by the chirp modulation as they 
have a wide support in signal space. The energy of these basis vectors is thus spread in fc-space thus improving the 
reconstruction quality of the low scale structures of the signal. On the other hand, the fine scale sparsity basis vectors 
at the center of the field of view remains not significantly modulated. However, these vectors are already incoherent 
with the Fourier basis as their energy naturally spreads out in fc-space. The fine scale structure of the signal are 
well recovered in absence and presence of chirp modulation. In summary, almost no shift-variant reconstruction 
quality is to be expected. 

Let us acknowledge that the idea of convolving the fc-space to optimize the acquisition procedure in the context 
of compressed sensing can also be found in [16]-[18]. In these works, the fc-space is convolved by a random 
Toeplitz matrices. We should also note that the spread spectrum technique can be related to the random convolution 
approach where the signal is convolved by a random sequence and under-sampled in real space [40]. In our case, 
convolution and under-sampling occur in k-space. Finally, let us mention that in a discrete setting, i.e., in the 
absence of band limit extension, replacing the linear chirp modulation by a random modulation leads to a universal 
encoding strategy, i.e., the reconstruction quality does not depend on the sparsity basis [20]. The S2MRI technique 
tries to emulate this universal encoding strategy. A universal compressed sensing strategy might as well be obtained 
by projecting the signal p onto random waveforms with Gaussian distributions [34], [35]. In the context of MRI, 
an encoding strategy based on this sensing scheme was recently proposed in [19]. 
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Fig. 3. Original brain image at 1 mm of resolution. The first panel represents its absolute value. The second panel represents its phase mapped 
between — tt and tt. The last two panels represent the logarithm of the amplitude of its Fourier transform before (third panel) and after (fourth 
panel) chirp modulation with chirp rate to = 0.3. Dark and light regions respectively indicate low and high intensities. 



III. Spread spectrum MRI (S2MRI) 
A. Quadratic phase profiles 

The spread spectrum principle was explained in the previous section in a simplified analog framework. Here, we 
move to a setting encompassing realistic analog MR signals which are limited in field of view and are consequently 
not band limited. 

Standard MR measurements take the form of fc-space coefficients, with k = {kx, ky, k^), of the original three- 
dimensional image probed representing the tissue magnetization as a function of the position inside a given field 
of view. This image is complex-valued due to magnetic field inhomogeneities [41]. We consequently denote it 
by a complex-valued function p{x) of the position x E R^, with components {x,y,z) in image space. (x,y) 
conventionally corresponds to the phase encoding directions, and z corresponds to the readout direction. The field 
of view is L — Lx X Ly X Lz. We also consider quadratic phase profiles represented by a linear chirp in the phase 
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encoding directions 

c(a;)=e''^K="'+"'''^'), (8) 

with chirp rate (wx,Wy) E M^. This linear chirp is characterized by an instantaneous frequency {wxX,Wyy) at 
position X. On the finite field of view L, it is therefore approximately a band-limited function with approximate 
band-limits {\wx\Lri./2, \wy\Ly/2) in the phase encoding directions. 

In this setting, MR measurements at spatial frequency fc e M"^ take the general form 

i/(fe)= / p{x)c{x)c-^''''''-'^d^x. (9) 

In other words, the measurement jy corresponds to the coefficient at spatial frequency fc of a signal obtained 
as the product of the original image p(x) with the linear chirp modulation c{x). In the absence of modulation 
(Wx — Wy ~ 0), the measurements simply reduce to standard fe-space measurements. In the presence of quadratic 
phase profiles, the modulation amounts to the convolution of the Fourier transform of the chirp with that of the 
original image. 

B. Under-sampling in k-space and the inverse problem 

Let us assume that we want to probe the signal p at a resolution corresponding to a band-limit B = {B^, By, Bz). 
Note that it does not imply that the signal is band-limited and some energy may remain beyond B. At this 
resolution, the signal of interest is discretized on a discrete uniform grid of A^ = N^ x Ny x A^^ spatial frequencies 
ki e M'', 1 < i < A^, in fc-space, with A^^; = '^L^B^, Ny = 2LyBy, and A^^ ~ 2LzBz- In real space, this 
discretized signal may equivalently be described by its coefficients on a discrete uniform grid of A^ points Xi 
with 1 < i < A^. On this discrete grid, the linear chirp can be parametrized in terms of discrete chirp rates 
{wx, Wy) ~ (wxL'^/Nx, WyLy/Ny'j, thus exhibiting approximate band-limits {\wx\Bx, \wy\By, 0) on the finite field 
of view L. 

We assume that the spatial frequencies probed span the fc-space up to the band-limits B. We also assume that 
these frequencies belong to the discrete grid of points hi, so that we can avoid any re-gridding operation [41]. Due 
to the linear chirp modulation, the measurements can contain significant energy from fc-space coefficients beyond 
the band-limits B. Considering the estimated band-limits of the linear chirp, we choose to reconstruct the original 
signal p on a high resolution grid of A^c = Nx (l + \wx\) x A^j, {l + \wy\)x Nz points in order to prevent any aliasing 
in the reconstruction algorithm from these high spatial frequencies. The sampled signal on this high resolution grid 
is denoted by a vector p E C^". The reconstructed signal is subsequently down-sampled at the desired resolution 
by keeping only its spatial frequencies belonging to the previously defined grid of size A^. 

The fc-space coverage provided by the M spatial frequencies probed feb, with 1 < 6 < M, can be represented by 
a binary mask in fc-space, equal to 1 for each spatial frequency probed and otherwise. The measurements may 
be denoted by a vector of M fc-space coefficients u = {z^(fc6)}i<h<M <= C^^ , possibly affected by noise, which 
is denoted by the vector n = {nb}i<6<M G C^ . In this setting, we consider an incomplete fc-space coverage 
(M < N), in order to accelerate the acquisition time in comparison with a complete fc-space coverage. In order to 
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Fig. 4. Simulated reconstructions for an under-sampling of 20 percent of A^ and an input snr of 32. The top row shows the reconstructions 
for the variable density sampling (left panel) and the S2MRI (right panel) techniques, respectively. The bottom row shows the error images 
(difference between the original image and the reconstructions) in the absence (left panel) and presence (right panel) of chirp modulation. For 
a better visualization, the error images were scaled by a factor of 6. The colormap for the error images goes from white to black, indicating 
low and high errors, respectively. 
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Fig. 5. Relative reconstruction errors as functions of the input snr for the variable density sampling technique (dashed blue curve) and the 
S2MRI technique with id = 0.3 (dot-dashed red curve) and w = iDopt (continuous black curve). The first to sixth panels show the curves for 
coverages of 5, 10, 15, 20, 25 and 50 percent of A^ respectively. All curves represent the mean relative error over 30 simulations, and the 
vertical lines represent the error at 1 standard deviation. 
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satisfy standard MRI constraints, an arbitrary under-sampling can be considered in the phase encoding directions 
while all spatial frequencies in the readout direction k^ are probed. This form of under-sampling can directly be 
expressed in terms of an acceleration of the acquisition. Thus, the S2MRI measurement model satisfies 

V ^MF*C\Jp + n, (10) 

where the matrix MF*CU G cMxAt^ gjjcodes the complete linear relation between the signal and the measurements. 
The rectangular matrix U G c^^uxAt^ represents an up-sampling operation, implemented in fc-space space by zero 
padding. This operation is needed to avoid any aliasing of the modulated signal due to a lack of sampling resolution 
in a discrete description of the originally continuous problem. The modulated signal is correctly obtained by applying 
the chirp modulation to the signal after up-sampling on the grid of N^ = Nx{l + 2|wa;|) x Ny{l + 2\'Wy\) x Nz 
points by zero padding in the k^ and ky directions. The matrix C G C^"^^" is the diagonal matrix implementing 
the chirp modulation on the up-sampled grid. The unitary matrix F G C^"^^" stands for the discrete Fourier basis 
on this high resolution grid. The matrix M G M^^^" is the rectangular binary matrix implementing the mask. It 
contains only one non-zero value on each line, at the index of the fc-space coefficient corresponding to each of the 
spatial frequencies probed. 

Regarding the reconstruction of the signal p, (10) represents the measurement constraint. We take a statistical point 
of view and consider independent Gaussian noise on each measurement. Considering a candidate reconstruction p, 
the residual noise reads as n = 1/ — MF*CUp. The noise level estimator, defined as twice the negative logarithm 
of the likelihood associated with p, reads as 

M |_ 
2/- N _ V^ \^b 



X 



™ I- 12 



6=1 



where the symbol | • | stands for the complex norm and cr^"'") for the standard deviation of each of the real and 
imaginary parts of the noise component rib. This noise level estimator follows a chi-square distribution with 2M 
degrees of freedom. Typically, this estimator should be minimized by a good candidate for the reconstruction. As 
explained in [33], the measurement constraint on the reconstruction may be defined as a bound x^(Pi i^) 1^ ^^ with 
e^ corresponding to some large percentile of the chi-square distribution. In the remainder of the paper, we choose 
the 99th percentile. 

C. Variable density sampling 

As discussed in [1], the method selecting the spatial frequencies fcb, with 1 < 6 < M, is critical to achieve 
good reconstruction. Their approach is based on the use of a variable density sampling method. The distributions, 
which are empirically identified to provide optimized reconstructions, are defined by a power law decaying function 
/(fc) == (1 — |fc|/|fcm|)P + (3 for some power p > and real constant /3 G [—1, 1] (|fep — k^ + ky and km is 
the highest spatial frequency of the fc-space domain to probe). In the present work, the spatial frequencies fc^ are 
selected using the method^ of [1]: the values /(fei) are first thresholded to restrict them to the interval [0, 1]; N 

'Toolbox available at http://www.stanford.edu/~mlustig/SparseMRI.html 
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independent binary random variables 6 (ki) taking value 1 with probability / (fe^) are then generated; and finally 
the mask M selecting the fc-space coefficients for which S (ki) = 1 is created. 

Let us highlight the importance of the value /3 on the actual shape of the variable density profile. For a fixed 
power p, this constant is computed beforehand in order to ensure that the number of measurements is, on average, 
equal to its target value M. Given a fixed number M of measurements, we denote pM the value of p for which 
/3 = 0. For a small value of p (p < pm), one would intuitively predict that the measurements spread all over the 
fc-space domain with the density of points higher at the center of fc-space. In fact, for p < pM, one has /3 < and 
a whole fc-space region at high frequency remains unprobed. On the contrary, we have /3 > for p > pM and the 
fc-space is probed with a non-zero probability at the edges. Consequently, p ~ pM is the first power that ensures 
that the entire fc-space domain is probed with a non-zero probability. Also note that for p > pM, the center of the 
fc-space is fully sampled, as /3 > 0, and that the size of the fully sampled region increases when p increases. 

In practice, the choice of the power p of the variable density profile is difficult and should be adapted with the 
number of measurements in order to obtain the best reconstruction qualities. In the absence of chirp modulation, we 
performed reconstructions for several values of p and noticed that p = pm leads to the best performance^. In the 
presence of a linear chirp modulation, the power spectrum of the original signal is spread and flattened. However, 
in the range of the chirp rates studied, the power spectrum of the modulated signal remains peaked at the origin 
(Fig. 3). We therefore also apply a variable density sampling and choose the power p to be pM- Empirically, this 
value also leads to the best reconstruction qualities. Note that choosing p — pm seems reasonable, as, given the 
spread of the information in the phase encoding directions, one wants to distribute the measurements over the entire 
sampling region of size B and also limit the size of the fully sampled region at the fc-space center. 

Let us acknowledge that recent theoretical results obtained in [42] support the choice of such profiles in MRI. 

D. Numerical simulations 

1) Simulation protocol: To test the proposed technique, a real brain image is acquired on a 7T short bore, actively 
shielded, MR scanner (Siemens, Erlangen, Germany). The subjects provided written informed consent prior to the 
imaging session, according to the guidelines of the local ethics committee. The parameters of the acquisition are 
as follows: L^x LyX L^ — 224 x 168 x 208 mm'' with a resolution of 0.5 x 0.5 x 4 mm''. The matrix size is thus 
N^xNyXN^ ^ 448 x 336 x 52. The standard clinical MPRAGE sequence is used with echo time TE = 3.53 ms, 
inversion time TI = 1.3 s, repetition time TR — 3.5 s, and bandwidth BW = 200 Hz. Note that, for the sake of 
simplicity, we restrict our analyses to one two-dimensional z-slice of the original three-dimensional acquisition with 
an under-sampling in both phase encoding directions k^ and ky. Also note that the image used is complex-valued. 

In order to model an analog acquisition scheme, the original image at a resolution of 0.5 mm is used to compute 
the measurements but the reconstruction is performed at a resolution of 1 mm. The reconstructed images are 
compared to the image obtained with a full acquisition at 1 mm of resolution (Fig. 3). 

*In practice, pm is determined by an iterative process. Starting from p = 0, we increase its value by 0.5 until we obtain /? > for the 
number of measurements considered. 
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The parameters of our analyses are as follows. Firstly, acquisitions are considered for various numbers M of 
complex measurements corresponding to coverages of 5, 10, 15, 20, 25 and 50 percent of N . Secondly, instrumental 
noise is also added to the measurements as independent and identically distributed zero-mean Gaussian noise. The 
corresponding standard deviation a is identical for all the frequencies probed and we consider values of input snr^ 
of 2^, with 1 < .7 < 6. Thirdly, the chirp modulation studied has the same chirp rate w in both phase encoding 
directions, i.e., w — Wx — Wy, with values in the range [0,0.5]. Fourthly, the signals are reconstructed by solving 
the BP problem where the Total Variation (TV) norm of the signal p is substituted for the ti norm**. This problem 
is solved thanks the Douglas Rachford algorithm [44], [45]. Note that the TV norm in combination with wavelet 
sparsity basis is, for example, used for reconstructing MRI signals from under-sampled fc-space in [1], [6], or [14]. 
For our problem, we tested several multiscale representations such as Daubechies wavelets, steerable wavelets, 
or curvelets, but the best reconstructions were obtained with the TV norm. Finally, for each value of M and 
input snr considered, 30 simulations are generated with independent noise and mask realizations, and the relative 
reconstruction errors \\p — p*||2/||p||2 are computed for the variable density sampling and S2MRI techniques. For 
each value of M and input snr, the discrete chirp rate w)opt that gives the smallest relative error on average over 
the 30 simulations is recorded^. 

2) Simulation results: The magnitudes of the reconstructed images obtained with the S2MRI and the variable 
density sampling techniques for an acceleration factor of 5 and an input snr of 32 are presented in Fig. 4, along 
with the corresponding error images (magnitudes of the complex-valued differences between the original image and 
reconstructed images). The relative errors of the reconstructions as functions of the input snr for the six coverages 
considered and for both methods are reported in Fig. 5. 

For acceleration factors larger than 2, the S2MRI technique with w — 0.3 provides better reconstruction than the 
variable density sampling technique with an improvement up to 0.05 of the relative error. Indeed, the relative error 
is, on average, lower in the presence of the chirp modulation. The corresponding standard deviations are also much 
smaller, indicating that the S2MRI technique is more stable'*'. At an acceleration factor of 2, the variable density 
sampling technique gives slightly better reconstructions than the S2MRI technique with w — 0.3. However, with 
w = Wopt = 0.1, the S2MRI technique provides relative errors similar to those obtained with the variable density 
sampling technique. These results suggest to reduce the chirp rate as the number of measurements increases. This 
is coherent with the fact that modulation is not needed in the limit of no under-sampling. 

'The snr is defined as the ratio between the mean value of the complex magnitude of the original signal and the standard deviation of the 
noise a. 

*The TV norm of a signal is defined as the (.1 norm of the magnitude of its gradient [32], [43]. Note that the recovery condition (7) does 
not hold with this norm. However, one can notice that the TV norm of a signal p is very similar to the l\ norm of its decomposition in the 
Haar wavelet basis. In the light of the preliminary results of Section II-C, one can thus hope to obtain an improvement of the reconstruction 
quality in presence of chirp modulation. 

'S2MRI toolbox available at http://lts2www.epfl.ch/people/gilles. 

'"Some leftover variability for the variable density sampling technique might still be removed by increasing the number of simulations. 
However, this would not modify the results and comparison with S2MRI. 
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When comparing the magnitudes of reconstructed images in Fig. 4, differences between both methods do not 
appear obvious. However, one can notice that fine details are recovered better with the S2MRI technique: the vessel 
(white spot) indicated by a blue arrow appears in the S2MRI reconstruction, but not in the variable density sampling 
reconstruction. The error images bring more information, and one can notice that the errors are smaller in the 
presence of chirp modulation. In particular, the low scale structures, rendered incoherent with the chirp modulation, 
are better recovered. 

IV. S2MRI IMPLEMENTATION AND EXPERIMENTAL RESULTS 
A. Implementation 

The S2MRI technique is tested on the 7T scanner described in Section Ill-D with 3D acquisitions of a phantom 
and a human brain. For the brain experiment, the subjects provided written informed consent prior to the imaging 
session, according to the guidelines of the local ethics committee. 

The chirp modulation is implemented through the use of a second order shim coil x^ — y^. In our implementation, 
the chirp rate varies linearly with the readout time t (or equivalently kz) and is proportional to the intensity of the 
quadratic magnetic field k: w{t) — Wx{t) = —'Wy{t) = ^nt/ir, where 7 is the gyromagnetic factor The maximum 
chirp rate w"^^^ is reached si t — TE + At/2, the mean chirp rate w™°^" ^[i — ^Yi, and the minimum chirp rate 
^min ^[i — 'pE _ At/2 (At is the readout duration). This chirp modulation can be introduced in the measurement 
model by modifying (10) as follows: 

i/ = MF*yCF*Up + n. (12) 

In the above equation, F* e C^u^^u ^j^^j p* 1^ ^WuXAfu implement the Fourier transform along the z-direction 
and the xy-directions, respectively. The matrix C e i£^n^-kn^ implements the chirp modulation. Signals are thus 
reconstructed on a grid of N^ = N^{1 + |wJP^^|) x Ny{l + |u)^^^|) x N^ points, and N^^ N^{1 + 2|w)^"^'=|) x 

This modulation can be decomposed as a quadratic phase modulation in the (x, y) planes with chirp rate w"^"''^'^, 
combined with a linear phase modulation in fc^. This linear phase modulation produces shifts of the original signal 
by an amount proportional to k(x^ ^ V^) along the z direction. The chirp modulation used is thus not ideal, as it 
creates distortions of the original object and complicates the measurement matrix. However, the energy of sparsity 
basis vectors is still spread by the main chirp modulation with chirp rate ■w^°^^\ and the previous conclusions 
based on theory and simulations should still hold. This will be confirmed by the numerical experiments of Section 
IV-B. Note also that the reconstructed images are free of any distortion as the complete effect of the modulation 
is modeled in (12). Let us remark that these distortions might be avoided with the use of RF pulses or dedicated 
coils applied only during phase encoding. 

As in Section III-D, the S2MRI technique is compared to the variable density sampling method with p = pM- 
Full acquisitions (M/N ~ 1) are performed both in the absence and in the presence of the chirp modulation. The 
number of phase encodings {kx,ky) is then reduced retrospectively by applying a mask on the complete data. 

February 25, 2013 DRAFT 




Fig. 6. Relative reconstruction errors as functions of the input snr for tlie variable density sampling technique (dashed blue curve) and the 
S2MRI technique (dot-dashed red curve) with varying chirp modulation. From left to right and top to bottom, the panels show the curves for 
coverages of 15, 25, 50 percent of A^ respectively. All curves represent the mean relative error over 5 simulations, and the vertical lines represent 
the error at 1 standard deviation. 

The noise level is evaluated directly on the under-sampled data available for reconstruction. For each {kx,ky) 
pair measured, all the frequencies k^ are probed so that the signal is available as a function of z. The level of the 
noise is estimated on probed pairs (/cj., ky) at positions z that do not contain any signal. This noise level is identical 
in the presence and absence of chirp modulation. 

B. Numerical validation 

In this section, we perform simulations using the acquisition scheme described above to confirm that, even though 
the implementation of the chirp modulation is not ideal, the reconstruction quality is still enhanced with the S2MRI 
technique. 

For this numerical experiment, a brain volume was acquired using the standard clinical MPRAGE sequence on 
a field of view of L^ — 243 mm, Ly — 176 mm and L^ = 256 mm, with a resolution of 1 mm in each direction 
(N^ = 243, Ny = 176, N^^ = 256). The echo time is TE = 4.59 ms, the inversion time TI = 1.5 s, the repetition 
time TR = 3.5 s, the bandwidth BW — 250 Hz. As in Section III-D, in order to model an analog acquisition 
scheme, the original image at a resolution of 1 mm is used to compute the measurements, but the reconstruction 
is performed at a resolution of 2 mm. The reconstructed images are compared to the image obtained with a full 
acquisition at 2 mm of resolution. 

The parameters of our experiment are as follows. Firstly, acquisitions are considered for various numbers M of 
complex measurements corresponding to coverages of 15, 25, 50 percent of N. Secondly, instrumental noise is also 
added to the measurements as independent and identically distributed zero-mean Gaussian noise. The corresponding 
standard deviation a is identical for all the frequencies probed and we consider values of input snr of 2^, with 

1 < J < 6. Thirdly, the simulated chirp modulation has a chirp rate varying linearly with k^: \wx\ G [0.12,0.30] 
and Iwyl e [0.09, 0.22]. These values for the discrete chirp rate correspond to those used during the real experiment 
performed hereafter at 1 mm of resolution. However, relatively to the band-limit at 2 mm of resolution, the spectrum 
is naturally more spread than relatively to the band-limit at 1mm of resolution. Therefore, for a reconstruction at 

2 mm of resolution, the spectrum does not need to be spread as much as for a reconstruction at 1 mm of resolution. 
We thus divided the values of the chirp rate by 2 in both dimensions. Fourthly, the signals are reconstructed by 

February 25, 2013 DRAFT 



16 



Variable density sampling 



S2MRI 





■ ^ 
i 










.<:j- 



'V. 









1 



^ .1 



vf ^ A^' 



Fig. 7. Simulated reconstructions with varying chirp modulation for an under-sampling of 50 percent of A'^ and an input siir of 32. The top row 
shows the reconstructions of an axial slice for the variable density sampling (left panel) and the S2MRI (right panel) techniques, respectively. 
The bottom row shows the error images (difference between the original image and the reconstructions) in the absence (left panel) and presence 
(right panel) of chirp modulation. For a better visualization, the error images were scaled by a factor of 8. The colormap for the error images 
goes from white to black, indicating low and high errors, respectively. 



solving the BP problem where the TV norm of the signal p is substituted for the ii norm. Finally, for each value 
of M and input snr considered, 5 simulations are generated with independent noise and mask realizations, and the 
relative reconstruction errors are computed for the variable density sampling and S2MRI techniques. 

The relative errors of the reconstructions as functions of the input snr for the three coverages considered and 
for both methods are reported in Fig. 6. The magnitudes of a reconstructed axial slice obtained with the S2MRI 
and the variable density sampling techniques for an acceleration factor of 5 and an input snr of 32 are presented 
in Fig. 7, along with the corresponding error images. Conclusions of Section III-D still hold with this acquisition 
scheme. The relative error of reconstruction is lower in the presence of the chirp for the three coverage considered. 
When comparing the error images in Fig. 7, one can once more notice that the errors are smaller with the S2MRI 
technique. 
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C. Experiments 

For the phantom experiment, a gradient echo sequence is used on a field of view of L^; = Lj, = L^ = 192 mm in 
the three directions with a resolution of 1 mm {N^ — Ny — N^ — 192). The echo time is TE = 6 ms, the repetition 
time TR = 10 ms, the bandwidth BW = 400 Hz, and the quadratic magnetic field intensity k = 3000 /iT/m^. The 
discrete chirp rates satisfy \wx\ = \wy\ G [0.23,0.36]. 

For the brain experiment, the standard clinical MPRAGE sequence is used on a field of view of L^ ~ 243 mm, 
Ly = 176 mm and L^ = 256 mm, with a resolution of 1 mm in each direction (Nx = 243, Ny ~ 176, N^ — 256). 
The echo time is TE = 4.59 ms, the inversion time TI = 1.5 s, the repetition time TR = 3.5 s, the bandwidth 
BW — 250 Hz, and the quadratic magnetic field intensity k ~ 4500ptT/m^. The discrete chirp rates satisfy 
\wx\ e [0.25,0.61] and \wy\ e [0.18,0.44]. 

Slices of the 3D reconstructions obtained with the S2MRI technique for 15, 25, 50 and 100 percent of phase 
encodings are presented in Fig. 8 for the phantom and brain acquisitions. The images obtained by inverse Fourier 
transform of the fully sampled fc-space (after correction of the distortions due to the small shifts in the z direction) 
are also presented. Fig. 9 provides a comparison of the reconstructed images obtained with the S2MRI and the 
variable density sampling techniques. In the aim of providing further insight on how each method preserves image 
resolution or, in other words, captures shape and magnitude of high frequency features, we also provide one- 
dimensional spatial profiles in Fig. 10. The white lines in Fig. 9 indicate the location of these profiles. 

Firstly, as one would expect, the visual reconstruction quality improves when the number of phase encodings M 
increases. In the limit of a coverage of 100 percent, we cannot identify any loss of details between the reconstructed 
images and the ones obtained by inverse Fourier transform. Moreover, the reconstructed image contains much less 
noise than the image obtained by inverse Fourier transform. Indeed, at a coverage of 100 percent, the problem (3) 
is essentially reduced to a denoising problem. 

Secondly, as in Section III-D, the differences between both methods do not appear at first glance when comparing 
only the magnitudes of reconstructed images in Fig. 9. Unfortunately here, the ground truth image is not accessible, 
so the corresponding error images cannot be displayed. However, a thorough visual inspection reveals that, for 
acceleration factors larger than 2, some fine details are better recovered by our approach. On the phantom recon- 
structions for acceleration factors of 6.7 and 4 (coverages of 15 and 25 percent, respectively), the separation between 
the two biggest circles is more visible. The shapes of the circles are also more curved. On the brain reconstructions 
for an acceleration factor of 6.7, the vessels at the center of image are still visible with our method but not with 
the variable density sampling technique. The cerebral cortex also appears sharper. For an acceleration factor of 2, 
the thin layer separating the two hemispheres of the brain remains more visible with the S2MRI technique. 

Thirdly, the reconstructed spatial profiles of the phantom presented in Fig. 10 show that the shape and magnitude 
of high frequency features are, more often, slightly better recovered with the S2MRI technique for acceleration 
factors larger than 4 (blue and red arrows indicate features better reconstructed with S2MRI and variable density 
sampling respectively). This improvement is much more significant on the brain data (see arrows), and holds for 
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j^gmg-y ^if^jiijjjYi ajj(j brain reconstructions from real experimental data with the S2MRI technique. The first to fourth columns showThe 
magnitude of the reconstructions from 15, 25, 50, and 100 percent of phase encodings respectively. The fifth column shows the reference 
images obtained by inverse Fourier transform (F.T.) of the fully sampled fc-space. 
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Fig. 9. Phantom and brain reconstructions from real experimental data with the S2MRI technique (first and third rows respectively) and 
the variable density sampling (VDS) technique (second and fourth rows respectively). The first to fourth columns show the magnitude of the 
reconstructions from 15, 25, 50, and 100 percent of phase encodings respectively. The fifth column shows the reference images obtained by 
inverse Fourier transform (F.T.) of the fully sampled fc-space. The white lines on these images indicate the location of the spatial profiles 
presented in Fig. 10 
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Fig. 10. Phantom and brain reconstruction profiles from real experimental data with the S2MRI technique (first and third rows respectively) and 
the variable density sampling (VDS) technique (second and fourth rows respectively). The first to fourth columns show the profile magnitude 
of the reconstructions (continuous red curve) from 15, 25, 50, and 100 percent of phase encodings respectively as well as the reference profiles 
(dotted black curve) obtained by inverse Fourier transform (FT.) of the fully sampled fc-space. 
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acceleration factors larger than 2. 

Finally, the distortions are correctly taken into account by the operator C, as both fully sampled slices of the 
phantom with and without chirp are identical (see Fig. 9). For the brain images, some differences remain due to 
small movements of the subject between both acquisitions. Let us remark that no fitting of the chirp parameters 
(center and rates) was performed to improve on the theoretical values. This highlights sufficient stability of the 
technique relative to parameter approximations. 

V. Conclusion and discussion 

We presented a spread spectrum technique (S2MRI) designed to accelerate MR acquisitions by compressed 
sensing. It consists of pre-modulating the original image by a linear chirp, which results from the application of 
quadratic phase profiles, and then performing random fc-space under-sampling. Non-linear algorithms promoting 
signal sparsity are then used for image reconstruction. 

In the context of compressed sensing theory, the effectiveness of the technique is supported by a decrease of 
coherence between the sensing and sparsity bases due to the pre -modulation. Simulations in a simplified analog 
setting confirm that the enhancement of the image reconstruction quality is linked to the evolution of the mutual 
coherence. The S2MRI technique was compared with the state-of-the-art variable density sampling using realistic 
numerical simulations and real acquisitions. Simulation results shows that the S2MRI technique performs slightly 
better than the variable density sampling technique in terms of relative reconstruction error for acceleration factors 
larger than 2. The chirp modulation was also implemented on a 7T scanner with the use of a second order 
shim coil. Simulations of this implementation confirms again the slight superiority of S2MRI. Visual inspection 
of reconstructions obtained from real acquisition of phantom and in vivo data also shows that this first (non-ideal) 
implementation provides slightly better reconstruction qualities than the variable density sampling method. The 
S2MRI technique thus outperforms the variable density sampling technique in terms of all the criteria used for 
evaluation. 

Regarding future evolutions of the S2MRI technique, an implementation of the linear chirp modulation with the 
use of RF pulses or dedicated coils could simplify the measurement scheme by applying the chirp modulation only 
during phase encoding in order to avoid object distortions, and, in turn, further enhance the reconstruction quality. 
Moreover, fitting the effective chirp center and rates on the basis of the data could improve the measurement model 
and result in better reconstructions. 

Let us also emphasize the potential interest of the S2MRI technique from an enhanced resolution perspective 
as, in the presence of the chirp modulation, the original image is reconstructed at a high resolution in order to 
avoid any aliasing problems. One can indeed consider reaching a higher spatial resolution for a fixed acquisition 
time without probing higher spatial frequencies in practice, which would require stronger gradient coils. In this 
context, any regularization approach adding a sparsity prior can help to synthesize spatial frequency information 
higher than that contained in the data. But the chirp modulation implies that high spatial frequency information is 
actually probed at lower frequencies. For illustration. Fig. 11 shows a reconstructed {x,y) slice for a coverage of 
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Fig. 11. Axial {x,y) reconstructed slice for a coverage of 100 percent in the presence of chirp modulation on a high resolution grid of size 
A^c = 393 X 255 (right panel) and after dowsampling on the original grid of size N = 243 X 176. 

100 percent on a high resolution grid A^c = 393 x 255 as well as the image obtained after downsampling on the 
grid A^ = 243 x 176. One can notice that the high resolution image provides sharper details with fewer aliasing 
artifacts. In particular, the vessels are better resolved. 
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